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Quantum Monte Carlo diagonalization for many-fermion systems 
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In this study we present an optimization method based on the quantum Monte Carlo diagonal- 
ization for many-fermion systems. Using the Hubbard-Stratonovich transformation, employed to 
decompose the interactions in terms of auxiliary fields, we expand the true ground-state wave func- 
tion. The ground-state wave function is written as a linear combination of the basis wave functions. 
The Hamiltonian is diagonalized to obtain the lowest energy state, using the variational principle 
within the selected subspace of the basis functions. This method is free from the difficulty known 
as the negative sign problem. We can optimize a wave function using two procedures. The first 
procedure is to increase the number of basis functions. The second improves each basis function 



through the operators, e 



using the Hubbard-Stratonovich decomposition. We present an al- 



gorithm for the Quantum Monte Carlo diagonalization method using a genetic algorithm and the 
renormalization method. We compute the ground-state energy and correlation functions of small 
clusters to compare with available data. 



PACS numbers: 74.20.-z, 71.10.Fd, 75.40.Mg 



I. INTRODUCTION 



The effect of the strong correlation between elec- 
trons is important for many quantum critical phenom- 
ena, such as unconventional superconductivity (SC) and 
the metal-insulator transition. Typical correlated elec- 
tron systems are high-temperature superconductors HI, 
[3, S 01, heavy fermions]!, 0, 0, 13] and organic 
conductors 0. Recently the mechanisms of superconduc- 
tivity in high-temperature superconductors and organic 
superconductors have been extensively studied using var- 
ious two-dimensional (2D) models of electronic interac- 
tions. Among them the 2D Hubbard model[l3| is the 
simplest and most fundamental model. This model has 
been studied intensively using numerical tools, such as 
the Quantum Monte Carlo method [H [11 [H, [11 [H, [H, 



,[li,[l,j2l|,j2 



Carlo method ;25','2l|27 
the two-leg ladder Hubbard model was also investi- 
gated with respect to the mechanism of high- temperature 
superconductivity Hi, [H, H [13, S [H, [S, [HI • 

The Quantum Monte Carlo (QMC) method is a numer- 
ical method employed to simulate the behavior of corre- 
lated electron systems. It is well known, however, that 
there are significant issues associated with the applica- 
tion to the QMC. First, the standard Metropolis (or heat 
bath) algorithm is associated with the negative sign prob- 
lem. Second, the convergence of the trial wave function 
is sometimes not monotonic, and further, is sometimes 
slow. In past studies workers have investigated the possi- 
bility of eliminating the negative sign problem 21, 22, 24i]. 
If the negative sign problem can be eliminated, the next 
task would be to improve the convergence of the simula- 



2p. 24 1, and the variational Monte 
U2i]2i,[33,[3l|,[32,[3ll. Recently, 



tion method. 

In this paper we present an optimization method based 
on Quantum Monte Carlo diagonalization (QMD or 
QMCD). The recent developments of high-performance 
computers have lead to the possibility of the simulation of 
correlated electron systems using diagonalization. Typi- 
cally, and as in this study, the ground-state wave function 
is defined as 



tH 



(1) 



where H is the Hamiltonian and ipQ is the initial one- 
particle state such as the Fermi sea. In the QMD method 
this wave function is written as a linear combination 
of the basis states, generated using the auxiliary field 
method based on the Hubbard-Stratonovich transforma- 
tion; that is 



(2) 



where (j)m are basis functions. In this work we have 
assumed a subspace with Nutates basis wave functions. 
From the variational principle, the coefficients {cm} are 
determined from the diagonalization of the Hamiltonian, 
to obtain the lowest energy state in the selected sub- 
space {0m}- Once the Cm coefficients are determined, 
the ground-state energy and other quantities are calcu- 
lated using this wave function. If the expectation values 
are not highly sensitive to the number of basis states, 
we can obtain the correct expectation values using an 
extrapolation in terms of the basis states at the limit 
Nstates oo. Howcver, a more reliable procedure must 
be employed when the change in the values at the limit 
is not monotonic. In this study results are compared to 



2 



results obtained from an exact diagonalization of small 
clusters, such as 4x4 and 6x2 lattices. 

In the following section, Section II, we briefly re- 
view the standard Quantum Monte Carlo simulation ap- 
proach. In Section III a discussion of the Quantum Monte 
Carlo diagonalization, and an extrapolation method to 
obtain the expectation values, are presented. Section 
IV is a discussion of the optimization procedure which 
employs the diagonalization method. All the results ob- 
tained in this study are compared to the exact and avail- 
able results of small systems in Section V. Finally, a sum- 
mary of the work presented in this paper is presented in 
Section VI. 



II. QUANTUM MONTE CARLO METHOD 

The method of Quantum Monte Carlo diagonalization 
lies in the QMC method. Thus it is appropriate to first 
outline the QMC method. The Hamiltonian is the Hub- 
bard model containing on-site Coulomb repulsion and is 
written as 

H = tij {cl„Cj„ + h.c.) + UY^ nj^Uji, (3) 

where c]^. {cja) is the creation (annihilation) operator of 

t 



an electron with spin a at the j-th site and rij^ 
tij is the transfer energy between the sites i and j . iij 
for the nearest-neighbor bonds. For all other cases tij = 
0. /7 is the on-site Coulomb energy. The number of sites 
is N and the linear dimension of the system is denoted 
as L. The energy unit is given by t and the number of 
electrons is denoted as A^e- 

In a Quantum Monte Carlo simulation, the ground 
state wave function is 



^ = e-"^Vo 



(4) 



where i/jq is the initial one-particle state represented by 
a Slater determinant. For large t, e~'^^ will project out 
the ground state from ipo. We write the Hamiltonian as 
H = K+V where K and V are the kinetic and interaction 
terms of the Hamiltonian in Eq.Q, respectively. The 
wave function in Eq.Q is written as 



-ArK„-ATV\M 



e 



)*Vo, (5) 



for T = At • M. Using the Hubbard-Stratonovich 
transformation [ill |43I . we have 



exp(-Art/7ii-fnix) = X! cxp(2asj(ni| - riij 



Si=±l 



(6) 



for (tanha)2 = tanh(ArC//4) or cosh(2a) = e'^^^/^ The 
wave function is expressed as a summation of the one- 
particle Slater determinants over all the configurations 



of the auxiliary fields Sj = ±1. The exponential operator 
is expressed as 

(7) 

where we have defined 

Bl{{s,{l)'}) = Q-^-K„^~v„{{s^m)^ (8) 

for 

= 2acr^Sj7ijo. - ij/Ar^n,,:,, (9) 



if„ = - ^ (4<xCjT + h.c). 
The ground-state wave function is 



(10) 



(11) 



where 0^ is a Slater determinant corresponding to a con- 
figuration TO = {si{i)} (i = 1, • • ■ , TV; € = 1, ■ • ■ , M) of 
the auxiliary fields: 



l[BUs,iM))---BUs^{mo 



(12) 



The coefficients Cm are constant real numbers: ci — C2 — 
■ ■ ■ . The initial state ■00 is a one-particle state. If elec- 
trons occupy the wave numbers fci, /C2, • • • , fcjv„ for each 

T I 

spin fj, V'o is given by the product V'oV'o where -i/jg is the 
matrix represented asp^ 



V 



\ 



J 



(13) 



is the number of electrons for spin a. In actual calcu- 
lations we can use a real representation where the matrix 
elements are cos{ki ■ rj) or sin(fci • rj). In the real-space 
representation, the matrix of Vo-({si}) is a diagonal ma- 
trix given as 

V;({sJ) = diag(2a(TSi - C/Ar/2, • • • , 2aasN - UAt/2). 

(14) 

The matrix elements of are 

{K^)ij — —t i,j are nearest neighbors 

= otherwise. (15) 
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(/)^ is an TV X No- matrix given by the product of the ma- 
trices e"^'^'^", e^" and ipQ. The inner product is thereby 
calculated as a determinant [2^. 

(CC>=det(Ct0-). (16) 

The expectation value of the quantity Q is evaluated as 



(Q) 



mn \ 



If Q is a bilinear operator for spin a, we have 



(17) 



)<,Odet(07/t^-'x) 



E™det(</>^V^)det((/.-'^t^,T'^) 
det(CV^)det(fc"t0-^) 

^ y , ,det((/)'"Uf )det 



(18) 



The expectation value with respect to the Slater deter- 
minants (0^(5o-(/2) is evaluated using the single-particle 
Green's function! la. |22| . 



(19) 



In the above expression, P^n = det((/)^(/)^)det(0^'^0~'^) 
can be regarded as the weighting factor to obtain the 
Monte Carlo samples. Since this quantity is not neces- 
sarily positive definite, the weighting factor should be 
l^rjinl; the resulting relationship is, 

{Qa) — ^ -^17171 {Q a) 7nn / ^ ^ Pmn 

mn mn 

= ^ \Pmn\sign{Pjnn){Q„),nn/ ^ \Pnin\sign{Pjnn) 



(20) 



where sign{a) = a/\a\ and 

{Qa^rnji — 



(21) 



This relation can be evaluated using a Monte Carlo pro- 
cedure if an appropriate algorithm, such as the Metropo- 
lis or heat bath method, is emplovedfi^. The summa- 
tion can be evaluated using appropriately defined Monte 
Carlo samples, 



^ ^ 7^T,mnS^9niPmn){Qa)r 

Tlh? Em„ sign{P^n) 



(22) 



where umc is the number of samples. The sign problem 
is an issue if the summation of sign{Pmn) vanishes within 
statistical errors. In this case it is indeed impossible to 
obtain definite expectation values. 



III. QUANTUM MONTE CARLO 
DIAGONALIZATION 



A. Diagonalization 

Quantum Monte Carlo diagonalization (QMD) is a 
method for the evaluation of (Qa) without the negative 
sign problem. The configuration space of the probabil- 
ity ll-Pmnll in Ea. (p2| is generally very strongly peaked. 
The sign problem lies in the distribution of Pmn in the 
configuration space. It is important to note that the 
distribution of the basis functions (to — 1, 2, • • • ) is 
uniform since Cm are constant numbers: ci = C2 = ■ ■ • . 
In the subspace {4>m}, selected from all configurations 
of auxiliary fields, the right-hand side of Eq. pT)) can be 
determined. However, the large number of basis states 
required to obtain accurate expectation values is beyond 
the current storage capacity of computers. Thus we use 
the variational principle to obtain the expectation values. 

From the variational principle. 



(Q) 



Emn '^rnC 



n \ V^m ^ 'i'n 



Emn ^mCn{4'm^n) 



(23) 



where Cm (to = 1, 2, • ■ • ) are variational parameters. In 
order to minimize the energy 



E = 



Emn CmCn{(l)niH(l)n) 
EmnCmC„((?im^n) 



(24) 



the equation dE/dcn = (n = 1, 2, • • ■ ) is solved for. 



m m 

If we set 



-m \ V^m 



)=0. 



{4>mH4>n), 



-^mn — {4^m4^n} ^ 



the eigen equation is 



Hu = EAu, 



(25) 



(26) 



(27) 



(28) 



for u = (ci, C2, ■ • ■ )*. Since 0m (w = 1, 2, ■ • ■ ) are not 
necessarily orthogonal, A is not a diagonal matrix. We 
diagonalize the Hamiltonian A~^H, and then calculate 
the expectation values of correlation functions with the 
ground state eigenvector; in general A^^H is not a sym- 
metric matrix. 

In order to optimize the wave function we must in- 
crease the number of basis states {(pm}- This can be 
simply accomplished through random sampling. For sys- 
tems of small sizes and small U , we can evaluate the 
expectation values from an extrapolation of the basis of 
randomly generated states. 
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B. Extrapolation 

In Quantum Monte Carlo simulations an extrapola- 
tion is performed to obtain the expectation values for 
the ground-state wave function. If M is large enough, the 
wave function in Eq. (|lip will approach the exact ground- 
state wave function, ipexactj as the number of basis func- 
tions, Nstates , is increased. If the number of basis func- 
tions is large enough, the wave function will approach, 
"fpexact, as M is increased. In either case the method em- 
ployed for the reliable extrapolation of the wave function 
is a key issue in calculating the expectation values. If the 
convergence is fast enough, the expectation values can be 
obtained from the extrapolation in terms of 1/ Nutates- 
Note that although the extrapolation in terms of 1/M, 
or the time step At, has often been employed in QMC 
calculations, however, a linear dependence for 1/M or 
At will not necessarily guarantee, an accurate extrapo- 
lated result. The variance method was recently proposed 
in variational and Quantum Monte Carlo simulations, 
where the extrapolation is performed as a function of the 
variance. An advantage of the variance method lies is 
that linearity is expected in some cases [U, |4^: 



{Q) - Qexact OC V, 

where v denotes the variance defined as 



(30) 



and Qexact is the expected exact value of the quantity Q. 

The following brief proof clearly shows that the energy 
in Eq. (|30p varies linearly. If we denote the exact ground- 
state wave function as and the excited states as V'i 
(j = 1, 2, • • • ), the wave function can be written as 



(31) 



where we assume that a and hi are real and satisfy - 
= 1. If it is assumed that Hipg = Egtpg and Hipi - 
Eitpi, the energy is found to be 

E = {H) 

= a^{'4>gHi;g) + 2a'^b,{iPiHil)g) + ^ bibjiiPiHi;, 



= a^Eg + Y^bibjii^iH^l^j) 

ij 

= a^Eg + Y,blE,. 



The deviation of E from Eg is 



(32) 



5E 



E-En 



{a'-\)Eg + Y.^fE, 

i 

b\{E,) - Eg) 



where 62 = 1.^2 ^^^^ ^ Ej^'^j/Ej^j- The 
variance w of 7? is also shown to be proportional to b^ if 
6^ is small. Since (H^) = a^Eg + b'^{Ef) where (Ef) = 
J2j b'jEj b'j, V is evaluated as 



SE 



SE 



+ •••}, (34) 



for a constant C. Hence if b is small it is found that 
5E 



Eg C 



(35) 



The other quantities can be found if Qg ~ {^gQipg), 
which leads to the result 

{Q)-Qg = -b'^Qg + 2a'^bi{'ip,QiPg) + ^b,bj{'ip,QiPj). 

i ij 

(36) 

If Q commutes with H, and tpi are eigenstates of Q, (Q) — 
Qg is proportional to 6^. 



{Q)-Qa = -b^iQa-{Q^)), 



(37) 



(29) ^here (Qi) = J2ibi{^iQ^i)/J2^bl, thus (Q) - Qg 



In the general case [H, Q] ^ 0, {Q) — Qg is not neces- 
sarily proportional to b^. However, if the matrix element 
(TpiQTpg) is negligible, we obtain 

{Q)-Qg = -b''Qg + J2b^bJ{^^Q^J) 

-b\Qg-y^bM^,Qi^i)iY.^i)-m 



'9 ^ 



(33) 



This shows that (Q) — Qg is proportional to the vari- 
ance V. Thus, if {ipiQipg) is small, we can perform an 
extrapolation using a linear fit to obtain the expectation 
values. We expect that this is the case for short-range 
correlation functions, since the local correlation may give 
rise to small effects in the orthogonality of ijji and i/'g, i-e. 
{tpiipg) = 0. Hence the evaluations of local quantities 
will be much easier than for the long-range correlation 
fimctions. 



IV. OPTIMIZATION IN QUANTUM MONTE 
CARLO DIAGONALIZATION 

A. Simplest algorithm 

The simplest procedure for optimizing the ground- 
state wave function is to increase the number of basis 
states {4'm} by random sampling. First, we set t and 
M, for example, t = 0.1, 0.2, • • • , and M = 20, 30, • • • . 
We denote the number of basis functions as N states ■ We 
start with Nstates = 100 ^ 300 and then increase up to 
2000 or 3000. This procedure can be outlined as follows: 
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Al. Generate the auxiliary fields Si {i = 1, • • ■ , iV) 
in B1{{si\)) randomly for ^ = 1,---,M for 

(to = I,-- - ,Nstates), and generate Nstates basis wave 
function {0m}. 

A2. Evaluate the matrices H„m — {4>mH(f)n) and 
Anm = {4>m4>n), and diagonalizc the matrix A^^H to 
obtain ip = 'Y^^Cm4>m- Then calculate the expectation 
values and the energy variance. 

A3. Repeat the procedure from Al after increasing the 
number of basis functions. 

For small systems this random method produces 
reliable energy results. The diagonalization plays an 
importance producing fast convergence. 

Failure of this simple method sometimes occurs as the 
system size is increased. The eigenfunction of A^^H can 
be localized when the off-diagonal elements are small, 
meaning that some components of are large and others 
are negligible. A quotient of localization in the configura- 
tion space can be defined. For example, the summation 
of \cmf except 0„ with large c„ is a candidate for such 
property, 

/ 

Qioc = ^|Cm|^, (39) 

m 

where the prime indicates that the summation is per- 
formed excluding the largest c„. Qioc should approach 1 
as the number of basis functions is increased. In the case 
of localization, Qioc < 0.1, where to lower the energy is 
procedurally inefficient. In order to avoid the localization 
difficulty there are two possible procedures. First is to 
multiply 4>rn. by B'^{{si})) to improve and optimize the 
basis wave function further. Second, use a more effec- 
tive method to generate new basis functions, explained 
further in the subsequent sections. 

B. Renormalization 

The basis functions {0™} multiplied by B'^ 
{£ — M + 1, M + 2, ■ ■ ■) are improved to provide a 
lower ground state. Here the 'improvement' means 
the increase of r in Eq.(|3]) which is accomplished by 
increasing M. The matrix B^{{si})) is given by a sum- 
mation over 2''^ configurations of {si}. If we consider 
all of these configurations, the space required for basis 
functions becomes large. Thus, we should select several 
configurations or one configuration that exhibits the 
lowest energy. One procedure to choose such a state is 
the following: 

Rl. Multiply (pm by exp(2acrsjnjcr — ^UArrija-), 
where we generate the auxiliary fields Sj {£) for ^ = M -I- 1 
and i = 1, ■ ■ ■ , N using random numbers. Then evaluate 
the ground state energy. If the energy is lower, is 
defined as a new and improved basis function. If we 
have a higher energy, (f>m remains unchanged. Repeat 



this procedure to lower the ground state energy twenty 
to fifty times. 

R2. Repeat above for m = 1, • • • , Nstates- 

R3. Multiply (f>m by the kinetic operator g-'^'^^r and 

R4. Repeat from Rl and continue tor £ ^ £ + 1. 

This method is referred to as the 1/2^-method in 
this paper since one configuration is chosen from 2^ 
possible states. It is important to note that Nstates 
remains unchanged. An alternative method has been 
proposed to renormalize {(j)m} and is outlined as (2^: 

R'l. Multiply (pm by Yia e:xp{2acrsjnjcr — ^UArrijcr) and 

evaluate the energy for sj = 1 and Sj = —1. We adopt 

Sj for which we have the lower energy. 

R'2. Repeat this procedure for j = I,-- - ,iV and 

determine the configuration {sj} for 0„i. 

R'3. Multiply by the kinetic operator e^^'^^^ and 

R'4. Repeat above for to ~ 1, • • • , Nstates to improve 
4>m, and repeat from Rl. 

In this latter method the energy is calculated for 
the auxiliary field s, = ±1 at each site before making a 
selection. In the literature [l^l this procedure is called 
the path-integral renormalization group (PIRG) method. 



C. Genetic algorithm 

In order to lower the ground-state energy efficiently, 
we can employ a genetic algorithm[4^ to generate the 
basis set from the initial basis set. One idea is to replace 
some parts of {si{£)} (i = 1, • ■ • , N;£ — 1, ■ ■ ■ , M) in 
that has the large weight |c„|^ to generate a new basis 
function 0^. The new basis function (j)'^ obtained in this 
way is expected to also have a large weight and contribute 
to ip. 

Let us consider two basis functions (pm and 4>n chosen 
from the basis set with a probability proportional 
to the weight \cj\ using uniform random numbers. 
For example, since J^aiij I'-il^ ~ 1' the weight 

of (pi to occupy I]j=i|cjl^ < a; < I]i=i|cjl^ the 
range < x < 1. If the random number r is within 

Sjl'iMcjl^ < r < X]j=i|cjl^; choose (pm, and cpn 
is similarly chosen. A certain part of the genetic data 
between (pm and (pn is exchanged, which results in two 
new basis functions cp'^ and (p'„. We add 0^, or (p'^, 
or both of them, to the set of basis functions as new 
elements. In this process every site is labeled using 
integers such as i = 1, • ■ • , iV, and then we exchange 
Si for i = Li,Li + 1, • • • , Li + L^xch — 1 where the 
number of Si to be exchanged is denoted as L^xch- Li 
can be determined using random numbers. We must 
also include a randomly generated new basis function as 
a mutation. Here we fix the numbers Nstates and Nstep 
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before starting the Monte Carlo steps. For instance, 
^states ^ 200 and N^tep = 200. Nutates IS increased 
as the Monte Carlo steps progress. We diagonalize the 
Hamiltonian A~^H at each step when the Nstep basis 
functions are added to the basis set in order to recal- 
culate the weight \ckf {k = 1, 2, ■ • • ). The procedure is 
summarized as follows: 

Gl. Generate the auxiliary fields Si{£) {i = 1, • • ■ ^N) 
randomly for = I,-- - .M. Generate N grates basis 
functions {(t)k\- This is the same as Al. 
G2. Evaluate the matrices Hmn = {4'mH(f>n) and 
^mn = i'Pm'Pn)) 8.nd diagonalizc the matrix A^^H to 
obtain tp = J2m ^^^m and calculate the expectation 
values and the energy variance. This is the same as A2. 
G3. Determine whether a new basis function should 
be generated randomly or using the genetic method on 
the basis of random numbers. Let Tc be in the range 

< < 1, for example, Vc = 0.9. If the random number 
r is less than rc, a new basis function is defined using 
the genetic algorithm and the next step G4 is executed, 
otherwise generate the auxiliary fields {sj} randomly 
and go to G6. 

G4. The weight of ipk is given as \ck\ ■ Choose two 
basis functions 0„i and from the basis set with 
a probability proportional to the weight \ckf. Now 
we determine which part of the genetic code is ex- 
changed between and We choose £ = £o for 

1 < £ < M using random numbers. We choose the sites 
j = Li, - ■ ■ , L2 = Li + Lexch — 1 for a randomly chosen 

G5. Exchange the genetic code {si{£)} between 0^ and 
(f>n for £ = £0 and j = Li, • • • , L2 + L^xch - 1- We have 
two new functions i/i^ and (j)'^. We adopt one or two of 
them as basis functions and keep the originals 0„i and 
(j)n in the basis set. 

G6. If the Nstep basis functions are added up to the basis 
set after step G2, then repeat from step G2, otherwise 
repeat from step G3. 



D. Hybrid optimization algorithm 

In actual calculations it is sometimes better to use 
a hybrid of genetic algorithm and renormalization 
method. The concept to reach the ground-state wave 
function employed in this study is presented in Fig.l. 
There are two possible paths; one is to increase the 
number of basis functions using the genetic algorithm 
and the other is to improve each basis function by the 
matrix Bi{{si}). The path followed when the hybrid 
procedure is employed is the average of these two 
paths and is represented as the diagonal illustrated in 
Fig.l. Before step G6 in the genetic algorithm, the 
basis functions (/),„ are multiplied by B£{{si}) following 
the renormalization algorithm of the steps Rl to R3. 



Optimization of each basis function 



Initial function 




FIG. 1: Concept of optimization procedure. There are three 
approaches to reach the ground-state wave function. First is 
to increase the number of basis functions for fixed m. Second 

is to increase Ad multiplying cacii basis function by Bi{{si}). 
Third is the hybrid method of the previous two procedures. 

Then we go to G6. The method is summarized as follows: 

HI. Generate the auxiliary fields Si{i) {i = I,-- - ,N) 
randomly for £ = .M. Generate N states basis 

hmctions {(pk}- 

H2. Evaluate the matrices Hmn = {(pmHcpn) and 
Amn = {<Pm(f>n), and diagonalize the matrix A~^H to 
obtain = ^^^^ Cm4>m and calculate the expectation 
values and the energy variance. 

H3. Determine whether a new basis should be generated 
randomly or using the genetic algorithm. Let rc be in 
the range < Tc < 1. If the random number r is less 
than rc, a new basis function is defined using the genetic 
algorithm and the next step is H4, otherwise generate 
the auxiliary fields {sj} randomly and go to H6. 
H4. The weight of (pk is given as |cfe|^. Choose two 
basis functions (f>m and 0„ from the basis set with 
a probability proportional to the weight |cfc|^. Now 
we determine which part of the genetic code is ex- 
changed between and We choose £ = £0 for 
1 < £ < M using random numbers. We choose the sites 
j = Li, - ■ ■ ,L2 = Li + Lexch — 1 for a randomly chosen 
Li. 

H5. Exchange the genetic code {si} between (p„i and 0„ 
for £ = £0 and j determined in step H4. We have two 
new functions (j)'^ and (j)'^. We adopt one or two of them 
as basis functions and keep the originals (pm and 0„ in 
the basis set. 

H6. Multiply 0m by exp(2acrsjnjcr — Aruja), 
where wc generate the auxiliary fields Si{£) for 
£ = M -t- land i = I,-- - ,N using random num- 
bers. Then evaluate the ground state energy. If the 
energy is lower, 4>„i is defined as a new and improved 
basis function. If we have a higher energy, 0^ remains 
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FIG. 2: Energy as a function of the variance for 4x 4, U = 4 
and A'^e = 10. The square is the exact result. The data 
fit using a straight line using the least-square method as the 
variance is reduced. We started with Nstates ~ 100 (first solid 
circle) and then increase up to 2000. 



FIG. 3: Energy as a function of the variance for 6 x 2 A^e = 10 
and U — 4. The square is the exact value obtained using exact 
diagonalization. 



unchanged. Repeat this procedure to lower the ground 

state energy twenty to fifty times. 

H7. Repeat above for m = 1, • • • , Nstates- 

H8. Multiply (pm by the kinetic operator g-^'^^r and 

H9. If the Nstep basis functions are added up to the 
basis set after step H2, then repeat from H2, otherwise 
repeat from step H3. 



E. Discussion on the Quantum Monte Carlo 
Diagonalization 

The purpose of the QMD method is to calculate 

In an algorithm based on the Quantum Monte Carlo pro- 
cedures, we evaluate the expectation values in the sub- 
space selected from all the configurations of the 
auxiliary fields. From the data showing how the mean 
values {Q) varies as the subspace is enlarged, we can es- 
timate the exact value of (Q) using an extrapolation. A 
devised algorithm may help us to perform the Quantum 
Monte Carlo evaluations efficiently. We have presented 
the genetic algorithm and the renormalization method. It 
may be possible to overcome the problem of localization 
in the subspace using this algorithm. In fact, the quo- 
tient Qioc in Eg. ([55)1 becomes nearly 1, i.e. Qioc > 0.99, 
in the evaluations presented in the next section. For such 
a case, most of basis functions in the subspace give con- 
tributions to the mean values of physical quantities and 
the obtained results are certainly reliable. 




-35 I ' ' ' ^ 

0.02 0.04 0.06 0.08 

Variance 



FIG. 4: Energy as a function of the variance v for 6x6. with 
the periodic boundary conditions. Solid circles and crosses 
are data obtained from the QMD method for two different 
initial configurations of the auxiliary fields. Gray open circles 
show results obtained from the 1/2'^-renormalization method 
(PIRG) with 300 basis wave functions. 



V. RESULTS 

In this section, the results obtained using the QMD 
method are compared to the exact and available results. 
We investigate the small clusters (such as 4 x 4 and 6 x 
6), the one-dimensional (ID) Hubbard model, the ladder 
Hubbard model, and the two-dimensional (2D) Hubbard 
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FIG. 5: Correlation functions obtained by QMD for 4 x 4 
lattice with Ne = 10 and U = 4 as a function of 1/Nstates- 



tions defined as follows. The Gutzwiller function is well 
known as 

V'G = Pg^O: (41) 
where Pq is the Gutzwiller projection operator, 

Pg = l[[l - {I -9)n,^n,^]. (42) 
j 

g is the parameter in the range < g < 1. The 
non-interacting wave function ipQ is optimized by con- 
trolling the double occupancy '^j{nj^nj\) . The fur- 
ther optimization of the Gutzwiller function can be 
obtainedfUSi, 



^(1) = e-^^e-"^ 



;(2) 



e-^'^e-"'^Vi'\ 



(43) 



(44) 



where K is the kinetic energy term and V is the on-site 
Coulomb interaction, 



model. 



A. Ground-state energy and correlation functions: 
check of the method 

The results for the 4x4, 6x2 and 6x6 systems are pre- 
sented in Table I. The results are compared to the exact 
values and those available values obtained using the ex- 
act diagonalization, the quantum Monte Carlo method, 
the constrained path Monte Carlo method [l^l and the 
variational Monte Carlo method for lattices with peri- 
odic boundary conditions. The expectation values for the 
ground state energy are presented for several values of U . 
The data include the cases for open shell structures where 
the highest-occupied energy levels are partially occupied 
by electrons. In the open shell cases the evaluations are 
sometimes extremely difficult. As is apparent from Table 
I, our method gives results in reasonable agreement with 
the exact values. The energy as a function of the variance 
is presented in Figs. 2, 3 and 4. To obtain these results 
the genetic algorithm was employed to produce the basis 
functions except the open symbols in Fig. 4. The 4x4 
where Ne = 10 in Fig. 2 is the energy for the closed shell 
case up to 2000 basis states. The other two figures are 
for open shell cases, where evaluations were performed up 
to 3000 states. Open symbols in Fig. 4 indicate the en- 
ergy obtained using the renormalization method (1/2''^- 
method) with 300 basis states. The results for the QMD 
and 1/2^-method (or PIRG) are quite similar as a func- 
tion of the energy variance. In these cases Qioc is close to 
1; Qioc ~ 0.99. As the variance is reduced, the data can 
fit using a straight line using the least-square method. 

In Table I we have also included the VMC results for 
the A-functions. The A-functions are variational func- 



V = 



(45) 



where A, a, A', a' are variational parameters to be deter- 
mined, to lower the ground-state energy, a is related to g 
as a = log(l/.g). This type of wave function is referred to 
as A- function in this paper. In our calculations the sec- 

(2) 

ond level A-function has given good results for the 
ground-state energy. If we perform an extrapolation as a 
function of the variance, we can obtain the correct expec- 
tation values as the QMD method. We must, however, 
determine variational parameters in the multi-parameter 
space by adjusting the values of the parameters to find a 
minimum. The advantage of the variational procedure is 
that the evaluations are stable even for large U/t, beyond 
the band width. 

The correlation functions for the 4x4 where — 10 
and U = 4 are presented in Table II. The exact diagonal- 
ization results are also provided. The correlation func- 
tions are defined as 



^(q) 



N ^ 



iq-(Rj-Ri) 



(("jT -"ji)('^»T -"ii)>: (46) 



s(i, j) = ((n^i - njx)(n,t - rii^)), (48) 

c{i,j)^{njni)-{nj){ni), (49) 

where rij = rij-^ + nji and Rj denotes the position of the 
j-th site. Aap is the pair correlation function, 



A„^(£) = (Ai(^-|-^)A^(^)), 



(50) 
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FIG. 6: Spin (solid circle) and charge (open circle) correla- 
tion functions obtained from the QMD method for the one- 
dimensional Hubbard model with 80 sites. The number of 
electrons is 66. We set U = 4 and use the periodic boundary 
condition. 



where Aa{i), a = x,y, denote the annihilation operators 
of the singlet electron pairs for the nearest-neighbor sites: 



(51) 



Here ci is a unit vector in the a(= a;, y)-direction. The 
agreement in this case is good for such a small system. 
The correlation functions are also dependent on the num- 
ber of basis wave functions as shown in FiglS] Since 
the fluctuation of the expectation values is small in this 
case, the extrapolation can be performed in terms of the 

ates ■ 



B. ID and Ladder Hubbard models 

In this subsection we show the results for the one- 
dimensional (ID) Hubbard model and ladder Hubbard 
model. The ground state of the ID Hubbard model is no 
longer Fermi liquid for U > 0. The ground state is insu- 
lating at half-filling and metallic for less than half-filling. 
The Fig. [S] is the spin and charge correlation functions, 
S{k) and C{k), as a function of the wave number, for 
the ID Hubbard model where N = 80. The 2kF sin- 
gularity can be clearly identified where the dotted line 
is for U = 0. The spin correlation is enhanced and the 
charge correlation function is suppressed slightly because 
of the Coulomb interaction. The momentum distribution 
function n(fc). 



n{k) 



5S( 



'fco- 



(52) 



FIG. 7: Momentum distribution function obtained from the 
QMD method for the one- dimensional Hubbard model with 
80 sites for the periodic boundary condition. The number of 
electrons is 66 and the Coulomb repulsion is U = 4. The 
dotted line is the guide given by ~ 0.5 -|- 0.4|fc — kF\^~^ 
where — 1 ~ 0.035 which corresponds to Kp ~ 0.69 using 
the formula 77 — 1 = {Kp + Kp ^ ) /4 — 1/2 [50l | . Open circles are 
the results obtained using the Gutzwiller function. 



is presented in FiglT] for the electron filling n = 0.825. 
Here Cka- is the Fourier transform of Cja- n{k) in the 
metallic phase exhibits a singular behavior near the wave 
number kp. The singularity close to fcj? is consistent with 
the property of the Luttinger hquidfsol [Slj . although it 
is difficult to analyze the singularity in more detail using 
the Monte Carlo method. The Gutzwiller function gives 
the unphysical result that n{k) increases as k approaches 
kp from above the Fermi surface. 

In the ladder Hubbard model, 



ladder 



= 1,2 ja 



JO- 



(53) 



where t{td) is the intrachain (interchain) transfer energy. 
The ladder Hubbard model exhibits a spin gap at half- 
filling, and the charge gap is also possibly opened for large 
[/ > at half-filling. The existence of superconducting 
phase has been suggested for the Hubbard ladder using 
the DMRG method 38] and the VMC method [11]. 

The spin correlation function ^(k) for the Hubbard 
ladder is presented in FigjSl where U = A and td = 1. 
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FIG. 8: Spin correlation function obtained from the QMD 
method for the ladder Hubbard model for 60 x 2 sites with 
periodic boundary condition. The number of electrons is 80 
and U = A. The upper line is for the upper band and the lower 
line is for the lower band. Singularities are at kpi — /cf2, 2fcF2, 
kpi + kF2 and 2A;_fi from left. The dotted lines are for U = Q. 
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FIG. 9: Momentum distribution function obtained from the 
QMD method for the ladder Hubbard model for 60 x 2 sites 
and periodic boundary condition. The number of electrons is 
80 and U = A. 



S{k) is defined as 

(54) 

where R^^ denotes the site (i,^) {£—1^2). We use the 
convention that k = (k,ky) where ky = and tt indicate 
the lower band and upper band, respectively. There are 



FIG. 10: Pair correlation function (solid circles) obtained 
using the QMD method for the ladder Hubbard model with 
16 X 2 sites where the boundary condition is open. (7 = 4, 
td = 1.4 and the electron filling is 0.875. The dashed line is 
the pair correlation function for [7 = 0. The open circles are 
the DMRG resuhs from Ref.fsl. 




FIG. 11: Charge gap as a function of U for = 1 (circles'). 
The DMRG results (squares) are provided for comparison [401] . 



four singularities at 2fci;'i, 2kF2, ^_fi — ^^2, and kFi+kF2 
for the Hubbard ladder, where kpi and kF2 are the Fermi 
wave numbers of the lower and upper band, respectively. 
They can be clearly identified as indicated by arrows in 
Figli 
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FIG. 12: Magnetization as a function of U for the half-filled 
Hubbard model after extrapolation at the limit of large A'^. 
Solid circles are the QMD results, and open circles are results 
obtained from the QMC methodT!^. The squares are the 
Gutzwiller-VMC results and gray solid circles show the 
3rd A- function {ijj'x^) VMC results carried out on the 8x8 
lattice[43l. The diamond symbol is the value from the two- 
dimensional Heisenberg model where m = 0.615[46l. l47j. 



The momentum distribution in Fig[5] 



(55) 



exhibits singularities at kpi and kp2 where the results 
obtained from the Gutzwiller function are also shown for 
comparison. Here we used the same notation for k and 
Ri^. The unphysical property of n(k) near the Fermi 
wave numbers for the Gutzwiller function are remedied 
in the QMD method. 

The pair correlation fimction, /S.yy{l) versus £ was also 
evaluated to compare with the density matrix renormal- 
ization group (DMRG) method. lS.yy{t) is defined as 



A,,(£) 



(56) 



FIG. 13: Momentum distribution function for the 14 x 14 
lattice. Parameters are i7 = 4 and A^e = 146. The boundary 
conditions are periodic in both directions. The results for the 
Gutzwiller function (open circle) are also provided. 



QMD method, we have presented the results for U = A. 
The enhancement of the pair correlation function over 
the non-interacting case is clear and is consistent with 
the DMRG method. 

It has been expected that the charge gap opens up as U 
turns on at half-filling for the Hubbard ladder model. In 
FigUDthe charge gap at half-filling is shown as a function 
of U . The charge gap is defined as 



Ae = E{N, + 2) + E{N, - 2) - 2E{N,), 



(58) 



for 



where E{Ne) is the ground state energy for the A^e elec- 
trons. The charge gap in Fig llll was estimated using the 
extrapolation to the infinite system from the data for the 
20 X 2, 30 X 2, and 40 x 2 systems. The data are consistent 
with the DMRG method and suggest the exponentially 
small charge gap for small U or the existence of the crit- 
ical value Uc in the range of < Uc < 1.5, below which 
the charge gap vanishes. 



2D Hubbard model 



Ay{i) = CiixC2iT - Cu^C2il. 



(57) 



Ayy{i) is the correlation function for the singlet pair on 
the rung. The results for Ayy {£) are given in FiglTO] on 
the 16 X 2 lattice for the open boundary condition, where 
the pair correlation functions Ayy (i) were averaged over 
several pairs, for a distance £. The values U — 4 and 
td — 1.4 are predefined, and the electron filling was n ~ 
0.875. The result obtained using the DMRG method is 
also provided for U = 8^3^ for comparison. Since a large 
value of U, such as J7 = 8, is not easily accessed using the 



The two-dimensional Hubbard model was also inves- 
tigated in this study. The results are presented in the 
following discussion. An important issue is the antifer- 
romagnetism at half-filling. The ground state is antifer- 
romagnetic for U > because of the nesting due to the 
commensurate vector Q = (tt, tt). The Gutzwiller func- 
tion predicts that the magnetization 



3 



hi. 



(59) 
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TABLE I: Ground state energy per site from the Hubbard 
model. The boundary conditions are periodic in both direc- 
tions. The current results are presented under the column 
labeled QMD. The constrained path Monte Carlo (CPMC) 
and Path integral renormalization group (PIRG) results are 
from Refs.Ull and H^], respectively. The column VMC is the 
results obtained for the optimized variational wave function 
t/)^^' except for 6x2 for which ipi^'' is employed. The QMC 
results are from Ref . . Exact results are obtained using 
diagonalization (isl] . 
Size Ne U QMD VMC CPMC PIRG QMC Exactto fc/2^ 



basis functions {(t>m} and diagonalize the Hamiltonian in 
this subspace. We can optimize the wave function by 
enlarging the subspace. The simplest way is to increase 
the number of basis functions by randomly generating 
auxiliary fields {si}. The wave function can be further 
improved by multiplying each (jjm by . Although the 
matrix in Eq.® generates 2^ new basis functions, we 
must select some states from them to keep the number of 
basis functions small. Within the subspace with the fixed 
number of basis functions, an extension of 1/2^-method 
method (fc = 1,2, 



■ ) is also possible. 



4x4 


10 


4 


-1.2237 


-1.221(1) 


-1.2238 


4x4 


14 


4 


-0.9836 


-0.977(1) 


-0.9831 


4x4 


14 


8 


-0.732(2) 


-0.727(1) 


-0.7281 


4x4 


14 


10 


-0.656(2) 


-0.650(1) 




4x4 


14 


12 


-0.610(4) 


-0.607(2) 


-0.606 


6x2 


10 


2 


-1.058(1) 


-1.040(1) 




6x2 


10 


4 


-0.873(1) 


-0.846(1) 




6x6 


34 


4 


-0.921(1) 


-0.910(2) 




6x6 


36 


4 


-0.859(2) 


-0.844(2) 





-0.920 
-0.8589 



-0.925 
■0.8608 



increases rapidly as U increases and approaches m — 1 
for large U. In Fig [T^] the QMD results are presented for 
TO as a function of U. The previous results obtained using 
the QMC method are plotted as open circles. The gray 
circles are for the A-function VMC method and squares 
are the Gutzwiller VMC data. Clearly, the magnetization 
is reduced considerably because of the fluctuations, and 
is smaller than the Gutzwiller VMC method by about 50 
percent. 

The Fig. [13] is the momentum distribution function 
n(k), 



-1.2238 We have proposed a genetic-algorithm based method 
-0.9840|;q generate the basis wave functions. The genetic algo- 
-u. 'troj.ji-j^j^ -^^jtjgiy used in solving problems to find the op- 
'll'g^ggtiiiiized solution in the space of large configuration num- 

. ' _„„J3ers. We make new basis functions from the functions 
-1.05807^ • 1 , . , . „ I i2 r 

_Q gyg-rwith large weightmg factors |c„| . New functions pro- 
duced in this way are expected to have large weighting 
factors. If the localization quotient Qioc in Eg. ([55)1 is 
not small, we can iterate the Monte Carlo steps without 
using the 1/2^-method. 

We have computed the energy and correlation func- 



(60) 



where the results for the Gutzwiller VMC and the QMD 
are indicated. The Gutzwiller function gives the results 
that n{k) increases as k approaches hp from above the 
Fermi surface. This is clearly unphysical. This flaw of 
the Gutzwiller function near the Fermi surface is not ob- 
served for the QMD result. 



VI. SUMMARY 

We have presented a Quantum Monte Carlo diagonal- 
ization method for a many-fermion system. We employ 
the Hubbard-Stratonovich transformation to decompose 
the interaction term as in the standard QMC method. 
We use this in an expansion of the true ground-state wave 
function. We have considered the truncated space of the 



tions for small lattices to compare with published data. 
The results obtained in this study are consistent with the 
published data. In the case of the open shell structures, 
evaluations are difficult in general and the convergence 
is not monotonic. In this case the subspace of the basis 
functions must be large to obtain the expectation values 
from the extrapolation procedure. 

As for the extrapolation, the expectation value (Q) 
may approach Qexact in a non-linear way, 



{N states) 



(61) 



for some exponent 9. We must evaluate 9 to obtain 
Qexact, from an extrapolation in terms of the iV^^ates- We 
may be able to use a derivative method where 9 is deter- 
mined so that the derivative d{Q) / dN states approaches 
as Nstates increases. In this paper we adopted the re- 
cently proposed energy- variance method [2^ llsj. For the 
energy and local quantities, we can expect {Q) — Qexact oc 
V for the variance v. It is expected that the long-range 
correlations are not trivial to calculate since the orthog- 
onality {ipiQipg) ~ should hold for the ground state ipg 
and excited states ipi. 
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